
#include "./mds_density.h"
#include <math.h>
#include <stdlib.h>

void mds_density_calcul( mds_frame *frame )
{
	
	const mds_vector    *box       = & frame->box;
	const mds_vector *   position  = frame->sets[MDS_POSITION];
	const size_t         n         = frame->num_particles;
	const size_t         nx        = 10;
	register size_t     i;
	int tab_dens[nx], ix;	
	Real dx;
	
	Real x_max=0, x_max2=0;
	
	x_max  = box->x;
	x_max2 = 0.5 * x_max;
	dx     = box->x/nx;
	
	/** erase file **/
	FILE *fichier2 = fopen( "density.txt", "wt" );
	if(fichier2)
	{
		fprintf( fichier2, "#\n" );
		fclose(fichier2);
	}
	
	/** initialize tab **/
	for( i=0; i < nx; i++ )
	{
		tab_dens[i] = 0;
		printf("%d",tab_dens[i]);
	}
	printf("\n");
	
	/** fill tab **/
	for( i=0; i < n; i++ )
	{
		
		const mds_vector *r_i = &position[i];
		const Real        x_i = r_i->x;
		
		/**test g(x)**/
		ix = (int) floor((x_i + x_max2) / dx);
		if( ix < 0 ) ix = 0;
		if(ix >= (int)nx) ix = nx-1;
		
		tab_dens[ix] = tab_dens[ix] + 1;
		
		printf("%d %d \n", ix, tab_dens[ix]);
		
		
		
	}
	
	fichier2 = fopen( "density.txt", "at" );
	if(fichier2)
	{
		for( i=0; i < nx; ++i )
		{
			fprintf( fichier2, "%g %d\n", -x_max2 + ( ((Real)i)+0.5) * dx, tab_dens[i] );
		}
		fclose(fichier2);
	}
	
}

/** function used by qsort **/
static 
int compare_Real(const void *pa, const void *pb)
{
	const Real *a_ptr = pa;
	const Real *b_ptr = pb;
	
	const Real a = *a_ptr;
	const Real b = *b_ptr;
	
	if(a<b){return -1;}
	else if(a>b){return 1;}
	else {return 0;}
	
}

void mds_nombre_accumulation( mds_frame *frame, Real sigma )
{
	const mds_vector    *box       = & frame->box;
	const mds_vector *   position  = frame->sets[MDS_POSITION];
	const size_t         n         = frame->num_particles;
	register size_t      i, j;
	int  nbre_voisin = 40, nbre_voisin2 = nbre_voisin / 2; /*warning nbre_voisin paire*/
	int  l = 0;
	
	Real  tab_position[n];
	Real  tab_accumulation[n+1];
	Real  der_f[n+1];
	Real  sigma_l = 1.3 * sigma; /*for each side*/
	Real  delta_x[nbre_voisin+1], delta_xmax[nbre_voisin+1];  /*distance between particle and neighbours and distance max*/
	Real  penteinf = 0, pentesup = 0;
	const Real  Lmin = -box->x * 0.5 + sigma_l;
	const Real  Lmax =  box->x * 0.5 - sigma_l;
	
	/**initialize delta_xmax**/
	for(i=0; i<=nbre_voisin; i++)
	{
		delta_xmax[i] = sigma_l;
	}
	
	/**create file and initialize**/
	FILE *fichier = fopen( "densityacc.txt", "wt" );
	if(fichier)
	{
		fprintf( fichier, "#\n" );
		fclose(fichier);
	}
	
	/** initialize tabs **/
	tab_accumulation[0] = 0;
	
	for( i=0; i < n; i++ )
	{	
		j = i + 1;
		
		const mds_vector *r_i = &position[i];
		
		tab_position[i] = r_i->x;
		tab_accumulation[j] = 0;
		der_f[i] = 0;
	}
	
	/**position array sorting in increasing order**/		
	qsort(tab_position, n, sizeof(Real), compare_Real);
	
	/**fill tabs**/
	for(i=0; i<n; i++)
	{
		j = i + 1;
		tab_accumulation[j] = tab_accumulation[j-1] + 1;
	}
	
	/**tab_voisin fill**/
	for(i=0; i<n; i++)
	{	
		if(tab_position[i]<Lmin || tab_position[i]>Lmax)
		{
			continue;
		}
		else
		{	
			/**voisins inférieur**/
			for(l=1; l<(nbre_voisin2+1); l++)
			{
				if(i<l){break;}
				else
				{	
					    delta_x[l] = tab_position[i] - tab_position[i-l];
						
						if(delta_x[l]<=delta_xmax[l])
						{
							penteinf += (tab_accumulation[i] - tab_accumulation[i-l]) / delta_x[l];
						}
				}
				/*printf("%g\n",delta_x[j]);*/
			}
			
			/**voisins supérieur**/
			for(l=1; l<(nbre_voisin2+1); l++)
			{
				if(i+l>=n)
				{
					break;
				}
				else
				{	
					delta_x[l] = tab_position[i+l] - tab_position[i];
					
					if(delta_x[l]<=delta_xmax[l])
					{
						pentesup += (tab_accumulation[i+l] - tab_accumulation[i]) / delta_x[l];
					}
				}
			}	
			penteinf = penteinf / nbre_voisin2;
			pentesup = pentesup / nbre_voisin2;
			
			der_f[i] = penteinf + pentesup;
		}
		/*printf("%g %g\n",tab_accumulation[i],tab_position[155]);*/
		
	}
	
	/**taylor**/
	/*for(i=0; i<n ;i++)
	{
		der_f[i] = (tab_accumulation[i+1] - tab_accumulation[i]) / (tab_position[i+1] - tab_position[i]);
	}*/
	
	
	/**print file**/
	fichier = fopen( "densityacc.txt", "at" );
	if(fichier)
	{
		for( i=0; i < n; i++ )
		{
			l = i + 1;
			fprintf( fichier, "%g %g %g\n", tab_position[i], tab_accumulation[l]/n, der_f[i] );
		}
		fclose(fichier);
	}
}

void mds_nombre_accumulation_radial( mds_frame *frame )
{
	const mds_vector *   position  = frame->sets[MDS_POSITION];
	const size_t         n         = frame->num_particles;
	register size_t      i, j;
	Real  tab_position[n];
	Real  tab_accumulation[n+1];

		
	/**create file and initialize**/
	FILE *fichier = fopen( "densityaccrad.txt", "wt" );
	if(fichier)
	{
		fprintf( fichier, "#\n" );
		fclose(fichier);
	}
	
	/** initialize tabs **/
	tab_accumulation[0] = 0;
	
	for( i=0; i < n; i++ )
	{	
		j = i + 1;
		
		const mds_vector *r_i = &position[i];
		const Real r = mds_vector_length( r_i );
		
		tab_position[i] = r;
		tab_accumulation[j] = 0;
	}
	
	/**position array sorting in increasing order**/		
	qsort(tab_position, n, sizeof(Real), compare_Real);
	
	/**fill tabs**/
	for(i=0; i<n; i++)
	{
		j = i + 1;
		tab_accumulation[j] = tab_accumulation[j-1] + 1;
	}
	
	/**print file**/
	fichier = fopen( "densityaccrad.txt", "at" );
	if(fichier)
	{
		for( i=0; i < n; i++ )
		{
			j = i + 1;
			fprintf( fichier, "%g %g\n", tab_position[i], tab_accumulation[j]/n);
		}
		fclose(fichier);
	}			
}

void mds_density_profil( mds_frame *frame )
{
	/*const mds_vector *   position  = frame->sets[MDS_POSITION];
	const size_t         n         = frame->num_particles;
	register size_t      i, j;
	Real  tab_position[n];
	Real  tab_accumulation[n+1];*/
	
		
	/**create file and initialize**/
	/*FILE *fichier = fopen( "densityaccrad.txt", "wt" );
	if(fichier)
	{
		fprintf( fichier, "#\n" );
		fclose(fichier);
	}	*/
}





